SUMMARY
This R code is used to estimate the
relationship between oxygen consumption (MO2) and ambient oxygen partial
pressure (PO2) in the Common galaxias (Galaxias maculatus). It
is also used to estimate the critical partial pressure of oxygen for
aerobic metabolism (Pcrit), which is commonly understood as the
threshold below which oxygen consumption rate can no longer be
sustained. The associated article is “The role of osmorespiratory
compromise in hypoxia tolerance of the purportedly oxyconforming teleost
Galaxias maculatus”.
AIM The article aims to test whether Galaxias maculatus can maintain oxygen consumption (MO2) as ambient PO2 falls, and if so, at what level it reaches critical partial pressure of oxygen for aerobic metabolism (Pcrit).
AUTHORS
To be added
AFFILIATIONS
To be added
AIM
To be added
These are the settings for the html output. We will use this to make out index file on Git
#kniter seetting
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE, # no warnings
cache = TRUE,# Cacheing to save time when kniting
tidy = TRUE
)
These are the R packages required for this script. You will need to install a package called pacman to run the p_load function.
# this installs and load packages
# need to install pacman
pacman::p_load("ggplot2",
"ggthemes",
"ggfortify",
"gtExtras",
"igraph",
"dagitty",
"ggdag",
"ggridges",
"gghalves",
"ggExtra",
"gridExtra",
"corrplot",
"RColorBrewer",
"gt",
"gtsummary",
"grid",
"plotly", # data visualisation
"tidyverse",
"janitor",
"readxl",
"broom",
"data.table",
"devtools",
"hms", # data tidy
"marginaleffects",
"brms",
"rstan",
"performance",
"emmeans",
"tidybayes",
"vegan",
"betareg",
"lme4",
"car",
"lmerTest",
"qqplotr",
"respirometry",
"mclust",
# modelling
"datawizard",
"SRS" # data manipulation
)
Here are some custom function used within this script.
calcSMR: authored by Chabot D. used to estimate SMR with several different methods Claireaux and Chabot (2016) DOI: doi:10.1111/jfb.12833
calcSMR = function(Y, q = c(0.1, 0.15, 0.2, 0.25, 0.3), G = 1:4) {
u = sort(Y)
the.Mclust <- Mclust(Y, G = G)
cl <- the.Mclust$classification
# sometimes, the class containing SMR is not called 1 the following
# presumes that when class 1 contains > 10% of cases, it contains SMR,
# otherwise we take class 2
cl2 <- as.data.frame(table(cl))
cl2$cl <- as.numeric(levels(cl2$cl))
valid <- cl2$Freq >= 0.1 * length(time)
the.cl <- min(cl2$cl[valid])
left.distr <- Y[the.Mclust$classification == the.cl]
mlnd = the.Mclust$parameters$mean[the.cl]
CVmlnd = sd(left.distr)/mlnd * 100
quant = quantile(Y, q)
low10 = mean(u[1:10])
low10pc = mean(u[6:(5 + round(0.1 * (length(u) - 5)))])
# remove 5 outliers, keep lowest 10% of the rest, average Herrmann & Enders
# 2000
return(list(mlnd = mlnd, quant = quant, low10 = low10, low10pc = low10pc, cl = cl,
CVmlnd = CVmlnd))
}
calcO2crit: authored by Chabot D. used to estimate O2crit (Pcript). Claireaux and Chabot (2016) DOI: doi:10.1111/jfb.12833
Note: O2 is assumed to be in percentage of dissolved oxygen (DO) to work
calcO2crit <- function(Data, SMR, lowestMO2 = NA, gapLimit = 4, max.nb.MO2.for.reg = 20) {
# AUTHOR: Denis Chabot, Institut Maurice-Lamontagne, DFO, Canada first
# version written in June 2009 last updated in January 2015
method = "LS_reg" # will become 'through_origin' if intercept is > 0
if (is.na(lowestMO2))
lowestMO2 = quantile(Data$MO2[Data$DO >= 80], p = 0.05)
# Step 1: identify points where MO2 is proportional to DO
geqSMR = Data$MO2 >= lowestMO2
pivotDO = min(Data$DO[geqSMR])
lethal = Data$DO < pivotDO
N_under_SMR = sum(lethal) # points available for regression?
final_N_under_SMR = lethal # some points may be removed at Step 4
lastMO2reg = Data$MO2[Data$DO == pivotDO] # last MO2 when regulating
if (N_under_SMR > 1)
theMod = lm(MO2 ~ DO, data = Data[lethal, ])
# Step 2, add one or more point at or above SMR 2A, when there are fewer
# than 3 valid points to calculate a regression
if (N_under_SMR < 3) {
missing = 3 - sum(lethal)
not.lethal = Data$DO[geqSMR]
DOlimit = max(sort(not.lethal)[1:missing]) # highest DO acceptable
# to reach a N of 3
addedPoints = Data$DO <= DOlimit
lethal = lethal | addedPoints
theMod = lm(MO2 ~ DO, data = Data[lethal, ])
}
# 2B, add pivotDO to the fit when Step 1 yielded 3 or more values?
if (N_under_SMR >= 3) {
lethalB = Data$DO <= pivotDO # has one more value than 'lethal'
regA = theMod
regB = lm(MO2 ~ DO, data = Data[lethalB, ])
large_slope_drop = (coef(regA)[2]/coef(regB)[2]) > 1.1 # arbitrary
large_DO_gap = (max(Data$DO[lethalB]) - max(Data$DO[lethal])) > gapLimit
tooSmallMO2 = lastMO2reg < SMR
if (!large_slope_drop & !large_DO_gap & !tooSmallMO2)
{
lethal = lethalB
theMod = regB
} # otherwise we do not accept the additional point
}
# Step 3 if the user wants to limit the number of points in the regression
if (!is.na(max.nb.MO2.for.reg) & sum(lethal) > max.nb.MO2.for.reg) {
Ranks = rank(Data$DO)
lethal = Ranks <= max.nb.MO2.for.reg
theMod = lm(MO2 ~ DO, data = Data[lethal, ])
final_N_under_SMR = max.nb.MO2.for.reg
}
# Step 4
predMO2 = as.numeric(predict(theMod, data.frame(DO = Data$DO)))
Data$delta = (Data$MO2 - predMO2)/predMO2 * 100 # residuals set to zero
# when below pivotDO
Data$delta[Data$DO < pivotDO | lethal] = 0
tol = 0 # any positive residual is unacceptable
HighValues = Data$delta > tol
Ranks = rank(-1 * Data$delta)
HighMO2 = HighValues & Ranks == min(Ranks) # keep largest residual
if (sum(HighValues) > 0)
{
nblethal = sum(lethal)
Data$W = NA
Data$W[lethal] = 1/nblethal
Data$W[HighMO2] = 1
theMod = lm(MO2 ~ DO, weight = W, data = Data[lethal | HighMO2, ])
# This new regression is always an improvement, but there can still
# be points above the line, so we repeat
predMO2_2 = as.numeric(predict(theMod, data.frame(DO = Data$DO)))
Data$delta2 = (Data$MO2 - predMO2_2)/predMO2_2 * 100
Data$delta2[Data$DO < pivotDO] = 0
tol = Data$delta2[HighMO2]
HighValues2 = Data$delta2 > tol
if (sum(HighValues2) > 0)
{
Ranks2 = rank(-1 * Data$delta2)
HighMO2_2 = HighValues2 & Ranks2 == 1 # keep the largest residual
nblethal = sum(lethal)
Data$W = NA
Data$W[lethal] = 1/nblethal
Data$W[HighMO2_2] = 1
theMod2 = lm(MO2 ~ DO, weight = W, data = Data[lethal | HighMO2_2,
])
# is new slope steeper than the old one?
if (theMod2$coef[2] > theMod$coef[2]) {
theMod = theMod2
HighMO2 = HighMO2_2
}
} # end second search for high value
} # end first search for high value
Coef = coefficients(theMod)
# Step 5, check for positive intercept
AboveOrigin = Coef[1] > 0
# if it is, we use a regression that goes through the origin
if (AboveOrigin) {
theMod = lm(MO2 ~ DO - 1, data = Data[lethal, ])
Coef = c(0, coefficients(theMod)) # need to add the intercept (0)
# manually to have a pair of coefficients
method = "through_origin"
HighMO2 = rep(FALSE, nrow(Data)) # did not use the additional value
# from Step 4
}
po2crit = as.numeric(round((SMR - Coef[1])/Coef[2], 1))
sum_mod = summary(theMod)
anov_mod = anova(theMod)
O2CRIT = list(o2crit = po2crit, SMR = SMR, Nb_MO2_conforming = N_under_SMR, Nb_MO2_conf_used = final_N_under_SMR,
High_MO2_required = sum(HighMO2) == 1, origData = Data, Method = method,
mod = theMod, r2 = sum_mod$r.squared, P = anov_mod$"Pr(>F)", lethalPoints = which(lethal),
AddedPoints = which(HighMO2))
} # end function
plotO2crit: used to plot the modes used for the calcO2crit function. Claireaux and Chabot (2016) DOI: doi:10.1111/jfb.12833
plotO2crit <- function(o2critobj, plotID = "", Xlab = "Dissolved oxygen (% sat.)",
Ylab = "dotitalumol", smr.cex = 0.9, o2crit.cex = 0.9, plotID.cex = 1.2, Transparency = T,
...) {
# AUTHOR: Denis Chabot, Institut Maurice-Lamontagne, DFO, Canada first
# version written in June 2009 last updated 2015-02-09 for R plotting
# devices that do not support transparency (e.g., postscript), set
# Transparency to FALSE
smr = o2critobj$SMR
if (Ylab %in% c("dotitalumol", "italumol", "dotumol", "umol", "dotitalmg", "italmg",
"dotmg", "mg")) {
switch(Ylab, dotitalumol = {
mo2.lab = expression(paste(italic(dot(M))[O[2]], " (", mu, "mol ", O[2],
" ", min^-1, " ", kg^-1, ")"))
}, italumol = {
mo2.lab = expression(paste(italic(M)[O[2]], " (", mu, "mol ", O[2], " ",
min^-1, " ", kg^-1, ")"))
}, dotumol = {
mo2.lab = expression(paste(dot(M)[O[2]], " (", mu, "mol ", O[2], " ",
min^-1, " ", kg^-1, ")"))
}, umol = {
mo2.lab = expression(paste(M[O[2]], " (", mu, "mol ", O[2], " ", min^-1,
" ", kg^-1, ")"))
}, dotitalmg = {
mo2.lab = expression(paste(italic(dot(M))[O[2]], " (mg ", O[2], " ",
h^-1, " ", kg^-1, ")"))
}, italmg = {
mo2.lab = expression(paste(italic(M)[O[2]], " (mg ", O[2], " ", h^-1,
" ", kg^-1, ")"))
}, dotmg = {
mo2.lab = expression(paste(dot(M)[O[2]], " (mg ", O[2], " ", h^-1, " ",
kg^-1, ")"))
}, mg = {
mo2.lab = expression(paste(M[O[2]], " (mg ", O[2], " ", h^-1, " ", kg^-1,
")"))
})
} else mo2.lab = Ylab
if (Transparency) {
Col = c(rgb(0, 0, 0, 0.7), "red", "orange")
} else {
Col = c(grey(0.3), "red", "orange")
}
Data = o2critobj$origData
lowestMO2 = quantile(Data$MO2[Data$DO >= 80], p = 0.05) # I added this
Data$Color = Col[1]
Data$Color[o2critobj$lethalPoints] = Col[2]
Data$Color[o2critobj$AddedPoints] = Col[3]
# ordinary LS regression without added points: blue line, red symbols
# ordinary LS regression with added points: blue line, red & orange symbols
# regression through origin: green dotted line, red symbols
line.color = ifelse(o2critobj$Method == "LS_reg", "blue", "darkgreen")
line.type = ifelse(o2critobj$Method == "LS_reg", 1, 3)
limX = c(0, max(Data$DO))
limY = c(0, max(Data$MO2))
plot(MO2 ~ DO, data = Data, xlim = limX, ylim = limY, col = Data$Color, xlab = Xlab,
ylab = mo2.lab, ...)
coord <- par("usr")
if (plotID != "") {
text(0, coord[4], plotID, cex = plotID.cex, adj = c(0, 1.2))
}
abline(h = lowestMO2, col = "pink") # I added this
abline(h = smr, col = "orange")
text(coord[1], smr, "SMR", adj = c(-0.1, 1.3), cex = smr.cex)
text(coord[1], smr, round(smr, 1), adj = c(-0.1, -0.3), cex = smr.cex)
if (!is.na(o2critobj$o2crit)) {
abline(o2critobj$mod, col = line.color, lty = line.type)
segments(o2critobj$o2crit, smr, o2critobj$o2crit, coord[3], col = line.color,
lwd = 1)
text(x = o2critobj$o2crit, y = 0, o2critobj$o2crit, col = line.color, cex = o2crit.cex,
adj = c(-0.1, 0.5))
}
} # end of function
meta_files_wd: Directory for the metadata
wd <- getwd()
meta_files_wd <- paste0(wd, "./meta-data") # creates a variable with the name of the wd we want to use
labchart_wd: Directory for Labchart estimated slopes
labchart_wd <- paste0(wd, "./lab-chart-slopes")
output_fig_wd: this is where we will put the figures
output_fig_wd <- paste0(wd, "./output-fig")
ifelse(!dir.exists("output-fig"), dir.create("output-fig"), "Folder already exists")
## [1] "Folder already exists"
labchart_df: We have imported the slopes extracted in LabChart during each phase of the experiment
setwd(labchart_wd)
#
# # Get the names of all sheets in the Excel file
sheet_names <- excel_sheets("labchart-all-dates_v2.xlsx")
all_trials_select <- c("start_date", "order", "phase", "cycle", "date", "time")
labchart_list <- list()
for (sheet in sheet_names) {
df <- read_excel("labchart-all-dates_v2.xlsx", sheet = sheet) %>%
dplyr::rename_with(tolower)
a_name <- paste0("a_", tolower(sheet))
a_df <- df %>%
dplyr::select(starts_with('a'), all_trials_select) %>%
dplyr::rename(temp = a_temp) %>%
dplyr::mutate(across(starts_with('a'), as.numeric)) %>%
pivot_longer(
cols = starts_with('a'), # Select all columns to pivot
names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
names_sep = "_"
) %>%
dplyr::mutate(respirometer_group = "a") # Add a new column with a fixed value
labchart_list[[a_name]]<- a_df
b_name <- paste0("b_", tolower(sheet))
b_df <- df %>%
dplyr::select(starts_with('b'), all_trials_select) %>%
dplyr::rename(temp = b_temp) %>%
dplyr::mutate(across(starts_with('b'), as.numeric)) %>%
pivot_longer(
cols = starts_with('b'), # Select all columns to pivot
names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
names_sep = "_"
) %>%
dplyr::mutate(respirometer_group = "b")
labchart_list[[b_name]] <- b_df
c_name <- paste0("c_", tolower(sheet))
c_df <- df %>%
dplyr::select(starts_with('c'), all_trials_select) %>%
dplyr::rename(temp = c_temp,
i_cycle = cycle) %>%
dplyr::mutate(across(starts_with('c'), as.numeric)) %>%
pivot_longer(
cols = starts_with('c'), # Select all columns to pivot
names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
names_sep = "_"
) %>%
dplyr::mutate(respirometer_group = "c") %>%
dplyr::rename(cycle = i_cycle)
labchart_list[[c_name]] <- c_df
d_name <- paste0("d_", tolower(sheet))
d_df <- df %>%
dplyr::select(starts_with('d'), all_trials_select) %>%
dplyr::rename(temp = d_temp,
i_date = date) %>%
dplyr::mutate(across(starts_with('d'), as.numeric)) %>%
pivot_longer(
cols = starts_with('d'), # Select all columns to pivot
names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
names_sep = "_"
) %>%
dplyr::mutate(respirometer_group = "d") %>%
dplyr::rename(date = i_date)
labchart_list[[d_name]] <- d_df
}
labchart_df <- bind_rows(labchart_list) %>%
dplyr::mutate(resp_cat_date = paste0(respirometer_group, "_", start_date),
chamber_n = str_extract(chamber_id, "\\d+"),
id_prox = paste0(resp_cat_date, "_", chamber_n),
time_hms = as_hms(time*3600),
date_chr = format(date, "%d/%m/%Y")
)
metadata: This is the meta data for each chamber
Note: We are also adding volume based on chamber type.
setwd(meta_files_wd)
metadata <- read_excel("Morpho.xlsx", na = "NA") %>%
dplyr::mutate(id_split = id) %>%
tidyr::separate(id_split, into = c("respirometer_group", "salinity_group", "start_date",
"chamber"), sep = "_") %>%
dplyr::mutate(volume = dplyr::case_when(chamber_type == "L" ~ 0.3, chamber_type ==
"M_M" ~ 0.105, chamber_type == "M_NM" ~ 0.11, chamber_type == "S" ~ 0.058,
chamber_type == "SM" ~ 0.075, chamber_type == "D3" ~ 0.055, TRUE ~ NA), id_prox = paste0(respirometer_group,
"_", start_date, "_", chamber))
Adding the meta data to LabChart slopes
labchart_tidy <- labchart_df %>%
dplyr::select(-start_date, -respirometer_group) %>%
left_join(metadata, by = "id_prox") %>%
dplyr::arrange(id)
We have 64 fish with MO2 data
n <- labchart_tidy %>%
dplyr::filter(chamber_condition == "fish") %>%
dplyr::distinct(id) %>%
nrow(.)
paste0("n = ", n)
## [1] "n = 64"
With 48 from the 0 ppt and 48 from 9 ppt groups
labchart_tidy %>%
dplyr::group_by(salinity_group) %>%
dplyr::reframe(`n total` = length(unique(id))) %>%
gt() %>%
cols_label(salinity_group = "Salinity group") %>%
cols_align(align = "center", columns = everything())
| Salinity group | n total |
|---|---|
| 0 | 48 |
| 9 | 48 |
Here we caculate the mean length and size of fish used in the experiment.
mass_length <- labchart_tidy %>%
dplyr::group_by(id) %>%
dplyr::sample_n(1) %>%
dplyr::ungroup() %>%
dplyr::reframe(x_mass = round(mean(mass), 3), min_mass = round(min(mass), 3),
max_mass = round(max(mass), 3), x_length = round(mean(length), 2), min_length = round(min(length),
2), max_length = round(max(length), 2))
mass_mean <- mass_length %>%
pull(x_mass)
mass_min <- mass_length %>%
pull(min_mass)
mass_max <- mass_length %>%
pull(max_mass)
length_mean <- mass_length %>%
pull(x_length)
length_min <- mass_length %>%
pull(min_length)
length_max <- mass_length %>%
pull(max_length)
paste0("The mean mass of fish was ", mass_mean, " g (range: ", mass_min, "–", mass_max,
")", ", and the mean length was ", length_mean, " mm (range: ", length_min, "–",
length_max, ")")
## [1] "The mean mass of fish was NA g (range: NA–NA), and the mean length was NA mm (range: NA–NA)"
We will remove 6 trials which had errors. These are as follows:
remove_trial_error <- c("a_0_25nov_3", "b_0_26nov_4", "c_0_22nov_2", "c_9_26nov_2",
"c_9_26nov_4", "d_9_27nov_3")
labchart_tidy <- labchart_tidy %>%
dplyr::filter(!(id %in% remove_trial_error))
We now have 58 fish with MO2 data
n <- labchart_tidy %>%
dplyr::filter(chamber_condition == "fish") %>%
dplyr::distinct(id) %>%
nrow(.)
paste0("n = ", n)
## [1] "n = 58"
With 45 in the 0 ppt group and 45 in the 9 ppt group
labchart_tidy %>%
dplyr::group_by(salinity_group) %>%
dplyr::reframe(`n total` = length(unique(id))) %>%
gt() %>%
cols_label(salinity_group = "Salinity group") %>%
cols_align(align = "center", columns = everything())
| Salinity group | n total |
|---|---|
| 0 | 45 |
| 9 | 45 |
Here we apply the following filters to the MO2 data:
cycle_burn <- 0:4
labchart_tidy <- labchart_tidy %>%
dplyr::filter(!(cycle %in% cycle_burn) & mo2corr < 0 & n > 60 & chamber_condition ==
"fish")
# Now we remove the points after the chamber is opened This is a function to do
# so
filter_o2_increase <- function(group) {
group <- group %>%
mutate(o2_diff = o2 - lag(o2)) # Calculate the difference in 'o2'
# Find the first index where 'o2_diff' exceeds 1
cutoff_index <- which(group$o2_diff > 1)[1]
# Filter the data up to the cutoff index, or return the full group if no
# cutoff
if (!is.na(cutoff_index)) {
group <- group[1:(cutoff_index - 1), ]
}
return(group)
}
# Apply the function to each group of 'chamber_id'
labchart_tidy_fish_closed <- labchart_tidy %>%
dplyr::filter(phase != "smr") %>%
group_by(id) %>%
group_split() %>%
lapply(filter_o2_increase) %>%
bind_rows() %>%
select(-o2_diff)
labchart_tidy_fish_smr <- labchart_tidy %>%
dplyr::filter(phase == "smr")
labchart_tidy_fish <- rbind(labchart_tidy_fish_smr, labchart_tidy_fish_closed) %>%
dplyr::arrange(id, order)
We have estimated SMR with two different appraches.
First using the mean of the lowest 3 values (smr_3l_means)
smr_3l_means <- labchart_tidy_fish %>%
dplyr::group_by(id) %>%
dplyr::filter(phase == "smr") %>%
dplyr::arrange(desc(mo2corr)) %>%
dplyr::slice_head(n = 3) %>% # Select the three lowest MO2
dplyr::ungroup() %>%
dplyr::group_by(id) %>%
dplyr::reframe(smr_l3 = mean(mo2corr))
# Combine the processed "smr" phase with all other phases
labchart_tidy_fish <- labchart_tidy_fish %>%
dplyr::left_join(., smr_3l_means, by = "id")
Second using the calcSMR function by Chabot, Steffensen and
Farrell (2016) DOI: 10.1111/jfb.12845. Specifically, We use mean of the
lowest normal distribution (MLND) where CVmlnd < 5.4 and the mean of
the lower 20% quantile (q0.2) were CVmlnd > 5.4. If CVmlnd is not
calculated we have used q0.2.
labchart_chabot_smr <- labchart_tidy_fish %>%
dplyr::filter(phase == "smr")
# Extract distinct IDs
ids <- labchart_chabot_smr %>%
dplyr::distinct(id) %>%
dplyr::pull()
# Initialise an empty list to store SMR data
smr_list <- list()
# Process each ID
for (id_i in ids) {
tryCatch({
# Filter the data for the current ID
df_i <- labchart_chabot_smr %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(abs_mo2corr = abs(mo2corr))
# Calculate SMR results
calcSMR_results <- calcSMR(df_i$abs_mo2corr)
CVmlnd_i <- calcSMR_results$CVmlnd
quant_i <- calcSMR_results$quant %>%
as_tibble()
quant_20per_i <- quant_i$value[3]
mlnd_i <- calcSMR_results$mlnd
smr_value <- if_else(CVmlnd_i < 5.4, mlnd_i, quant_20per_i)
smr_type <- if_else(CVmlnd_i < 5.4, "mlnd", "quant_20per")
smr_value <- if_else(is.na(smr_value), quant_20per_i, smr_value)
smr_type <- if_else(is.na(smr_type), "quant_20per", smr_type)
# Create a data frame for the current ID
smr_df <- tibble(id = id_i, smr = smr_value, smr_est = smr_type)
}, error = function(e) {
# Handle errors by assigning NA values
smr_df <- tibble(id = id_i, smr = NA, smr_est = NA)
})
# Append to the list
smr_list[[id_i]] <- smr_df
}
# Combine all individual SMR data frames into one
smr_df <- bind_rows(smr_list) %>%
dplyr::rename(smr_chabot = smr, smr_chabot_method = smr_est)
labchart_tidy_fish <- labchart_tidy_fish %>%
dplyr::left_join(., smr_df, by = "id")
Here we are transforming the MO2 units. The resulting vaules are as follows:
# Combine back into one data frame
labchart_tidy_fish <- labchart_tidy_fish %>%
dplyr::mutate(DO = conv_o2(
o2 = o2,
from = "mg_per_l",
to = "percent_a.s.",
temp = temp, #C
sal = measured_salinity,
atm_pres = 1013.25),
net_volume = volume - mass/1000,
MO2 = abs(mo2corr)*net_volume*60*60,
MO2_g = MO2/mass,
SMR = abs(smr_l3)*net_volume*60*60,
SMR_g = SMR/mass,
SMR_CHABOT = abs(smr_chabot)*net_volume*60*60,
SMR_CHABOT_g = SMR_CHABOT/mass
)
Here we plot all oxygen consumption (MO2; mg O2/g/h) by dissolved oxygen percentage (DO) for all fish, including all SMR estimates.
labchart_tidy_fish %>%
ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
geom_point(show.legend = FALSE) +
geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
geom_smooth(se = TRUE, colour = "red", linetype = "dashed") +
theme_clean() +
labs(
subtitle = "All values",
x = "Dissolved oxygen percentage (DO)",
y = "MO2 (mg O2 g/h)"
)
Same plot but without SMR values.
labchart_tidy_fish %>%
dplyr::filter(phase != "smr") %>%
ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
geom_point(show.legend = FALSE) +
geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
geom_smooth(se = TRUE, colour = "red", linetype = "dashed") +
theme_clean() +
labs(
subtitle = "Only closed periods",
x = "Dissolved oxygen percentage (DO)",
y = "MO2 (O2 mg/g/h)"
)
Looking at the difference responses in the two salinity groups.
It’s appears more variable in freshwater.
labchart_tidy_fish %>%
ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
geom_point(show.legend = FALSE) +
geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
geom_smooth(se = TRUE, colour = "red", linetype = "dashed") +
theme_clean() +
facet_wrap(~salinity_group) +
labs(
subtitle = "mo2 vs o2 by salinity treatment",
x = "Dissolved oxygen percentage (DO)",
y = "MO2 (O2 mg/g/h)"
)
Looking at the difference chamber types
labchart_tidy_fish %>%
ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
geom_point(show.legend = FALSE) +
geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
geom_smooth(se = TRUE, colour = "red", linetype = "dashed") +
theme_clean() +
facet_wrap(~chamber_type, scale = "free") +
labs(
subtitle = "mo2 vs o2 by chamber type",
x = "Dissolved oxygen percentage (DO)",
y = "MO2 (O2 mg/g/h)"
)
min_o2 <- min(labchart_tidy_fish$o2, na.rm = TRUE)
max_o2 <- max(labchart_tidy_fish$o2, na.rm = TRUE)
labchart_tidy_fish <- labchart_tidy_fish %>%
mutate(o2_group = cut(o2,
breaks = seq(min_o2, max_o2, length.out = 11), # 12 intervals, so 13 breakpoints
labels = paste0("Group ", 1:10),
include.lowest = TRUE))
time_bin_df <- labchart_tidy_fish %>%
dplyr::group_by(o2_group) %>%
dplyr::reframe(mean_MO2_g = mean(MO2_g)*31.25,
mean_o2 = mean(o2),
n = length(MO2_g)*31.25,
MO2_g_sd = sd(MO2_g)*31.25,
o2_sd = sd(o2))
time_bin_df %>%
ggplot(aes(y = mean_MO2_g, x = mean_o2)) +
# Add raw data points
geom_point(data = labchart_tidy_fish, aes(y = MO2_g*31.25, x = o2),
size = 2, color = "grey", alpha = 0.5) + # Raw data points
# Add summary points
geom_point(size = 3, colour = "black", show.legend = FALSE) +
# Add vertical error bars
geom_errorbar(aes(ymin = mean_MO2_g - MO2_g_sd, ymax = mean_MO2_g + MO2_g_sd),
width = 0.15, colour = "black") +
# Add horizontal error bars
geom_errorbarh(aes(xmin = mean_o2 - o2_sd, xmax = mean_o2 + o2_sd),
height = 0.4, colour = "black") +
theme_clean() +
labs(
subtitle = "All values with error bars",
x = "O2 (mg/L)",
y = "MO2 (umol O2 g/h)"
)
Plotting MO2 estimates for each fish. The dashed red line is
Chabot SMR methods, and the solid line is the mean of the lowest 3
measures (excluding the first 5 cycles)
Notes: There’s something wired going on with a_0_25nov_2 it seems like many of the raw MO2 values are positive.
# Create output directory if needed
output_fig_slopes_wd <- file.path(output_fig_wd, "slopes")
if (!dir.exists(output_fig_slopes_wd)) {
dir.create(output_fig_slopes_wd)
}
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
pull(id) %>%
as.list()
MO2_plot_list <- list()
# 1) Open the PDF device once
pdf(file = file.path(output_fig_slopes_wd, "combined_slopes.pdf"), width = 8, height = 6)
# 2) Loop over IDs and create each plot
for (id_i in ids) {
smr_chabot <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
dplyr::slice(1) %>%
dplyr::pull(SMR_CHABOT)
smr_l3 <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
dplyr::slice(1) %>%
dplyr::pull(SMR)
plot <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
ggplot(aes(x = o2, y = MO2)) + geom_hline(yintercept = smr_chabot, linetype = "dashed",
color = "darkred") + geom_hline(yintercept = smr_l3, color = "darkred") +
geom_point(aes(colour = phase)) + theme_clean() + labs(subtitle = paste0(id_i,
" slopes"), x = "Mean o2 (mg_per_l)", y = "abs(mo2) (mg_per_l)")
# Instead of saving each plot separately, just print it
print(plot)
MO2_plot_list[[id_i]] <- plot
}
# 3) Close the PDF device *after* the loop
dev.off()
## png
## 2
for (p in MO2_plot_list) {
print(p)
}
Here we are following the methods Urbina et al. (2012) with an
incremental regression analyses, in order to determine the best fit for
the data.
This analysis evaluated each polynomial order equation starting at zero and then increasing to the third order. This permitted a mathematical assessment of whether the data best fitted a single linear relationship (0th-order polynomial; suggesting the fish were oxyconforming and do not reach a Pcrit), or whether a PO2 crit value could be determined as an intersection point of two distinct functions (one at hypoxic oxygen concentrations, the other at normoxic; i.e. oxyregulation).
We will use weighted regression to account for a higher density of data at normoxic conditions (i.e. SMR values). Specifically, we are weighting the importance of each data point in fitting the model based on frequency of points in that space. Points that have higher weights influence the model fit more, while points with lower weights have less impact. A high density of points at high o2 values could lead to overfitting in that region, while underfitting or misrepresenting trends in lower-density regions (e.g., low o2 vaules). The weighting is achieved by dividing the o2 values into bins, computing the frequency of points in each bin, and assigning weights as the inverse of frequency for each point.
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
pull(id) %>%
as.list()
weighted_model_comparison_list <- list()
weighted_model_results_list <- list()
for (id_i in ids) {
# Filter data for the current ID
df_i <- labchart_tidy_fish %>%
dplyr::filter(id == id_i)
# Calculate weights based on O2 density
df_i <- df_i %>%
dplyr::mutate(
o2_bin = cut(o2, breaks = 10) # Divide O2 into 10 bins
) %>%
dplyr::group_by(o2_bin) %>%
dplyr::mutate(
bin_freq = length(order), # Count points in each bin
weight = 1 / bin_freq # Weight = inverse frequency
) %>%
ungroup()
# Fit models with weights
models <- list(
lm_0 = lm(MO2 ~ 1, data = df_i, weights = weight), # 0th-order (constant mean)
lm_1 = lm(MO2 ~ o2, data = df_i, weights = weight), # 1st-order (linear)
lm_2 = lm(MO2 ~ poly(o2, 2), data = df_i, weights = weight), # 2nd-order (quadratic)
lm_3 = lm(MO2 ~ poly(o2, 3), data = df_i, weights = weight) # 3rd-order (cubic)
)
# Extract metrics to compare models
weighted_model_comparison_list[[id_i]] <- purrr::map_df(models, glance, .id = "model") %>%
dplyr::mutate(id = id_i) %>%
dplyr::select(id, everything())
weighted_model_results_list[[id_i]] <- map_df(models, ~ tidy(.x, conf.int = TRUE), .id = "model") %>%
clean_names() %>%
mutate(id = id_i) %>%
select(id, everything())
}
# Combine dataframes into a single data frame
weighted_model_comparison <- bind_rows(weighted_model_comparison_list) %>%
dplyr::mutate(poly = as.numeric(stringr::str_remove_all(model, "lm_")))
weighted_model_results <- bind_rows(weighted_model_results_list) %>%
dplyr::mutate(poly = as.numeric(stringr::str_remove_all(model, "lm_")))
Selecting the best fitting model
best_weighted_model <- weighted_model_comparison %>%
dplyr::group_by(id) %>%
dplyr::arrange(desc(AIC)) %>%
dplyr::slice(1) %>%
dplyr::ungroup()
Most often the best fitting model is a 0th-order polynomial
(n = 32, 55.17%), suggesting that MO2 does not show a
statistically significant dependence on the o2. In other words, the
metabolic rate does not adjust based on oxygen availability, and there
is no clear critical oxygen threshold (Pcrit) where the relationship
changes. This is indicative of a oxyregulator.
The next most common is a 1st-order polynomial (n = 12, 20.69%) which may suggest the presences of linear relationship between o2 and MO2 in some fish, however, to be true evidence of a oxyconformer this relationship should be positive (i.e. as o2 falls MO2 also falls). Only 6 of the 12 individuals best modelled with a linear function were positive, and none were statistically significant (Table 1).
The reaming fish (n = 14, 24.14%) had a 2nd- or 3rd-oder polynomial which could suggest the presence of a critical oxygen threshold (Pcrit) where the relationship between o2 and MO2 changes. These would also need to be validated to check the same of the polynomials to deterimine if a Pcirt is present. B
total_fish <- nrow(best_weighted_model)
best_weighted_model %>%
dplyr::group_by(poly) %>%
dplyr::reframe(n = length(id), percent = round((n/total_fish) * 100, 2)) %>%
gt() %>%
cols_align(align = "center", columns = everything())
| poly | n | percent |
|---|---|---|
| 0 | 32 | 55.17 |
| 1 | 12 | 20.69 |
| 2 | 4 | 6.90 |
| 3 | 10 | 17.24 |
For the 12 fish that had o2 and MO2 relationships best modelled
with a linear function, only 6 are were positive relationships (which we
would expect if fish were oxyconforming), and none are statistically
significant.
ids_lm_1_list <- best_weighted_model %>%
dplyr::filter(model == "lm_1") %>%
dplyr::pull(id)
weighted_model_results %>%
dplyr::filter(id %in% ids_lm_1_list & model == "lm_1" & term == "o2") %>%
dplyr::mutate(estimate = round(estimate, 4), ci = paste0(round(conf_low, 4),
"–", round(conf_high, 4)), p_value = round(p_value, 3)) %>%
dplyr::select(id, estimate, ci, p_value) %>%
gt() %>%
cols_align(align = "center", columns = everything())
| id | estimate | ci | p_value |
|---|---|---|---|
| a_0_24nov_1 | -0.0006 | -0.0033–0.0021 | 0.648 |
| a_0_24nov_3 | 0.0017 | -0.0032–0.0066 | 0.495 |
| a_0_26nov_1 | -0.0009 | -0.0048–0.0029 | 0.618 |
| a_0_27nov_4 | -0.0009 | -0.0033–0.0014 | 0.435 |
| b_0_26nov_2 | 0.0044 | -0.0027–0.0115 | 0.214 |
| b_9_21nov_4 | 0.0009 | -0.0016–0.0034 | 0.471 |
| b_9_22nov_2 | 0.0002 | -9e-04–0.0013 | 0.728 |
| b_9_22nov_4 | -0.0037 | -0.0125–0.0051 | 0.400 |
| c_9_24nov_4 | -0.0002 | -9e-04–5e-04 | 0.551 |
| d_0_22nov_2 | 0.0009 | -9e-04–0.0027 | 0.297 |
| d_9_24nov_2 | -0.0012 | -0.0031–8e-04 | 0.233 |
| d_9_26nov_3 | 0.0001 | -0.001–0.0011 | 0.907 |
Visualising the regressions
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
pull(id) %>%
as.list()
for (id_i in ids) {
df_i <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(o2_bin = cut(o2, breaks = 10)) %>%
dplyr::group_by(o2_bin) %>%
dplyr::mutate(
bin_freq = length(order), # Count points in each bin
weight = 1 / bin_freq # Weight = inverse frequency
) %>%
ungroup()
best_weighted_model_i <- best_weighted_model %>%
dplyr::filter(id == id_i)
poly_i <- best_weighted_model_i %>%
dplyr::pull(poly)
poly_i_name <- best_weighted_model_i %>%
dplyr::mutate(name = case_when(
poly == 0 ~ "0th-order polynomial",
poly == 1 ~ "1st-order polynomial",
poly == 2 ~ "2nd-order polynomial",
poly == 3 ~ "3rd-order polynomial",
TRUE ~ "ERROR"
)) %>%
dplyr::pull(name)
r <- best_weighted_model_i %>%
dplyr::pull(r.squared) %>%
round(., 2)
sigma <- best_weighted_model_i %>%
dplyr::pull(sigma) %>%
round(., 2)
mean_MO2 <- df_i %>%
dplyr::reframe(mean = mean(MO2), na.rm = TRUE) %>%
dplyr::pull(mean)
x_max <- df_i %>%
dplyr::reframe(max = max(o2), na.rm = TRUE) %>%
dplyr::pull(max)
y_max <- df_i %>%
dplyr::reframe(max = max(MO2), na.rm = TRUE) %>%
dplyr::pull(max)
plot <- df_i %>%
ggplot(aes(x = o2, y = MO2, weight = weight)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ poly(x, poly_i), se = FALSE, colour = "blue") +
geom_hline(yintercept = mean_MO2, colour = "grey", linetype = "dashed", linewidth = 1) +
annotate("text", x = x_max/2,
y = y_max,
label = paste0(poly_i_name, "\n", "R = ", r, " Sigma = ", sigma),
hjust = 0, vjust = 1, size = 4) +
labs(
title = paste0(id_i),
x = "O2 (mg/L)",
y = "MO2 (mg O2 g/h)"
) +
theme_minimal()
print(plot)
}
Here we will calculate Pcrit using Chabot method and function
calcO2crit. We are using our estimates for SMR (mean of lowest three).
This function uses the fifth percentile of the MO2 values observed at
dissolved oxygen levels ≥ 80% air saturation as the criterion to assess
low MO2 values. The algorithm then identifies all the MO2 measurements
greater than this minimally acceptable MO2 value. Within this sub-set,
it identifies the ̇ MO2 measurement made at the lowest DO and thereafter
considers this DO as candidate for breakpoint (named pivotDO in the
script). A regression is then calculated using observations at DO levels
< pivotDO, and a first estimate of O2crit is calculated as the
intersection of this regression line with the horizontal line
representing SMR. The script then goes through validation steps to
ensure that the slope of the regression is not so low that the line,
projected to normoxic DO levels, passes below any MO2 values observed in
normoxia. It also ensures that the intercept is not greater than zero.
Corrective measures are taken if such problems are encountered.
lowestMO2 default is the quantile(Data\(MO2[Data\)DO >= 80], p=0.05). It is used to segment the data and locate the pivotDO.
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
dplyr::pull()
pcrit_model_df_list <- list()
pcrit_models <- list()
for (id_i in ids) {
df_i <- labchart_tidy_fish %>%
dplyr::filter(id == id_i)
o2crit <- calcO2crit(Data = df_i, SMR = df_i$SMR[1], lowestMO2 = NA, gapLimit = 4,
max.nb.MO2.for.reg = 7)
vaule <- o2crit$o2crit
lowestMO2 = quantile(df_i$MO2[df_i$DO >= 80], p = 0.05)
SMR <- o2crit$SMR
nb_mo2_conforming <- o2crit$Nb_MO2_conforming
r2 <- o2crit$r2
method <- o2crit$Method
p <- o2crit$P[1]
pcrit_model_df <- tibble(id = id_i, pcrit_vaule = vaule, pcrit_smr = SMR, pcrit_lowestMO2 = lowestMO2,
pcrit_nb_mo2_conforming = nb_mo2_conforming, pcrit_r2 = r2, pcrit_method = method,
pcrit_p = p)
pcrit_model_df_list[[id_i]] <- pcrit_model_df
pcrit_models[[id_i]] <- o2crit
}
pcrit_model_df <- bind_rows(pcrit_model_df_list)
Here’s the plots for the Pcrit estimates
# Create output directory if needed
output_fig_pcrit_chabot_wd <- file.path(output_fig_wd, "model_chabot")
if (!dir.exists(output_fig_pcrit_chabot_wd)) {
dir.create(output_fig_pcrit_chabot_wd)
}
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
dplyr::pull()
pcrit_chabot_list <- list()
# Open a single PDF device
pdf(file = file.path(output_fig_pcrit_chabot_wd, "combined_chabot_plots.pdf"), width = 8,
height = 6)
for (id_i in ids) {
r2 <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
dplyr::pull(pcrit_r2)
conforming <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
dplyr::pull(pcrit_nb_mo2_conforming)
P <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
dplyr::pull(pcrit_p)
SMR <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
dplyr::pull(pcrit_smr)
lowestMO2 <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
dplyr::pull(pcrit_lowestMO2)
# Generate and render the plot
plotO2crit(o2critobj = pcrit_models[[id_i]])
# Add a title
mtext(text = paste0(id_i), side = 3, line = 2, adj = 0, col = "blue", font = 2,
cex = 1.2)
mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
font = 1, cex = 0.8)
}
# Close the PDF device *after* the loop
dev.off()
## png
## 2
Printing in htlm document
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
dplyr::pull()
for (id_i in ids) {
comment <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
dplyr::slice(1) %>%
dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
pull(comment)
r2 <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
dplyr::pull(pcrit_r2)
conforming <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
dplyr::pull(pcrit_nb_mo2_conforming)
P <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
dplyr::pull(pcrit_p)
SMR <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
dplyr::pull(pcrit_smr)
lowestMO2 <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
dplyr::pull(pcrit_lowestMO2)
# Generate and render the plot
plotO2crit(o2critobj = pcrit_models[[id_i]])
# Add a title
mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
font = 2, cex = 1.2)
mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
font = 1, cex = 0.8)
}
We need to set some rules as to when the Pcrit estimates are
reliable, as it seems many of our fish do not seem to reach a Pcrit.
We can filter for only cases were at the lowest O2 value three consecutive MO2 measures full below our SMR and fifth percentile of the MO2 values observed at dissolved O2 levels >80%. In the model output these are called nb_mo2_conforming points.
pcrit_list <- pcrit_model_df %>%
dplyr::filter(pcrit_nb_mo2_conforming > 2) %>%
pull(id)
mean_pcrit <- pcrit_model_df %>%
dplyr::filter(pcrit_nb_mo2_conforming > 2) %>%
dplyr::reframe(mean = mean(pcrit_vaule)) %>%
pull(mean)
paste0("There are ", length(pcrit_list), " fish with identified Pcrits. ", "The mean Pcrit is ",
round(mean_pcrit, 2))
## [1] "There are 14 fish with identified Pcrits. The mean Pcrit is 22.83"
for (id_i in pcrit_list) {
comment <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
dplyr::slice(1) %>%
dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
pull(comment)
r2 <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
dplyr::pull(pcrit_r2)
conforming <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
dplyr::pull(pcrit_nb_mo2_conforming)
P <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
dplyr::pull(pcrit_p)
SMR <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
dplyr::pull(pcrit_smr)
lowestMO2 <- pcrit_model_df %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
dplyr::pull(pcrit_lowestMO2)
# Generate and render the plot
plotO2crit(o2critobj = pcrit_models[[id_i]])
# Add a title
mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
font = 2, cex = 1.2)
mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
font = 1, cex = 0.8)
}
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
dplyr::pull()
pcrit_model_df_list_2 <- list()
pcrit_models_2 <- list()
for (id_i in ids) {
df_i <- labchart_tidy_fish %>%
dplyr::filter(id == id_i)
o2crit <- calcO2crit(Data = df_i, SMR = df_i$SMR_CHABOT[1], lowestMO2 = NA, gapLimit = 4,
max.nb.MO2.for.reg = 7)
lowestMO2 = quantile(df_i$MO2[df_i$DO >= 80], p = 0.05)
vaule <- o2crit$o2crit
SMR <- o2crit$SMR
nb_mo2_conforming <- o2crit$Nb_MO2_conforming
r2 <- o2crit$r2
method <- o2crit$Method
p <- o2crit$P[1]
pcrit_model_df <- tibble(id = id_i, pcrit_vaule = vaule, pcrit_SMR = SMR, pcrit_lowestMO2 = lowestMO2,
pcrit_nb_mo2_conforming = nb_mo2_conforming, pcrit_r2 = r2, pcrit_method = method,
pcrit_p = p)
pcrit_model_df_list_2[[id_i]] <- pcrit_model_df
pcrit_models_2[[id_i]] <- o2crit
}
pcrit_model_df_2 <- bind_rows(pcrit_model_df_list_2)
Plotting with the SMR Chabot method
ids <- labchart_tidy_fish %>%
dplyr::distinct(id) %>%
dplyr::pull()
for (id_i in ids) {
comment <- labchart_tidy_fish %>%
dplyr::filter(id == id_i) %>%
dplyr::slice(1) %>%
dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
pull(comment)
r2 <- pcrit_model_df_2 %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
dplyr::pull(pcrit_r2)
conforming <- pcrit_model_df_2 %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
dplyr::pull(pcrit_nb_mo2_conforming)
P <- pcrit_model_df_2 %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
dplyr::pull(pcrit_p)
SMR <- pcrit_model_df_2 %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_SMR = round(pcrit_SMR, 3)) %>%
dplyr::pull(pcrit_SMR)
lowestMO2 <- pcrit_model_df_2 %>%
dplyr::filter(id == id_i) %>%
dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
dplyr::pull(pcrit_lowestMO2)
# Generate and render the plot
plotO2crit(o2critobj = pcrit_models_2[[id_i]])
# Add a title
mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
font = 2, cex = 1.2)
mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
font = 1, cex = 0.8)
}
Here using the 100 closed trials we will estimate Pcrit with
five popular techniques for Pcrit calculation: the traditional
breakpoint metric (broken stick regression), the nonlinear regression
metric (Marshall et al. 2013), the sub-prediction interval metric (Birk
et al. 2019), the alpha-based Pcrit method (Seibel et al. 2021), and the
linear low O2 (LLO) method (Reemeyer & Rees 2019).
The function is called calc_pcrit() and is part of the respirometry
package.
Link: https://search.r-project.org/CRAN/refmans/respirometry/html/calc_pcrit.html
Parameters to consider
avg_top_n: for alpha method, a numeric value
representing the number of top α0 (MO2/PO2) values to average together
to estimate α. Default is 1. We recommend no more than 3 to avoid
diminishing the α value with sub-maximal observations.
level: for Sub_PI method, Percentage at which
the prediction interval should be constructed.
iqr: Only for Sub_PI. Removes mo2 observations
that are this many interquartile ranges away from the mean value for the
oxyregulating portion of the trial. If this filtering is not desired,
set to infinity.
NLR_m: only applies to NLR. Pcrit is defined as
the PO2 at which the slope of the best fitting function equals NLR_m
(after the MO2 data are normalized to the 90% quantile). Default is
0.065
MR: A numeric value for the metabolic rate at
which pcrit_alpha and pcrit_LLO should be returned. If not supplied by
the user, then the mean MO2 of the “oxyregulating” portion of the curve
is applied for pcrit_alpha and NA is returned for pcrit_LLO.
mo2_threshold: A single numeric value above which mo2 values are ignored for alpha Pcrit estimation. Useful to removing obviously erroneous values. Default is Inf.
We will only use 100 c trails for this.
# labchart_tidy_fish_100c <- labchart_tidy_fish %>% dplyr::filter(phase ==
# '100c') labchart_tidy_fish_100c_n <- labchart_tidy_fish_100c %>%
# dplyr::distinct(id) %>% nrow(.) paste0('n for 100 closed = ',
# labchart_tidy_fish_100c_n)
Here we build the models
# combined_pcirt_list <- list() ids <- labchart_tidy_fish_100c %>%
# dplyr::distinct(id) %>% pull(id) %>% as.list() for (id_i in ids) { id_name <-
# id_i mo2_data <- labchart_tidy_fish_100c %>% dplyr::filter(id == id_i) MR_set
# <- mo2_data$SMR[1] %>% as.numeric() # Use tryCatch to handle errors and skip
# problematic calculations pcrit_df <- tryCatch({ pcrit_df <- calc_pcrit(po2 =
# mo2_data$o2, mo2 = mo2_data$MO2, method = 'All', avg_top_n = 2, # alpha
# metric (default = 1) recommend no more than 3 level = 0.95, # Sub_PI metric
# (default = 0.95) iqr = 1.5, # Sub_PI metric (default = 1.5) NLR_m = 0.065, #
# NLR metric (default = 0.065) MR = MR_set, # alpha and LLO metrics,
# mo2_threshold = Inf, # alpha metric return_models = FALSE # return model
# parameters? ) %>% as.data.frame() %>% rownames_as_column(var = 'method') %>%
# rename(value = '.') %>% tidyr::pivot_wider(., names_from = method,
# values_from = value) %>% dplyr::mutate(id = id_name) %>% dplyr::select(id,
# everything()) }, error = function(e) { message('Skipping channel ', id_name,
# ' due to error: ', conditionMessage(e)) NULL }) # Only add to list if
# pcrit_df is not NULL if (!is.null(pcrit_df)) { combined_pcirt_list[[id_name]]
# <- pcrit_df } }
Combined all the Pcrit model estimates together
# pcirt <- bind_rows(combined_pcirt_list)
Here we will save the plots for the various Pcrit curves.
# # Create output directory if needed output_fig_pcrit_100c_wd <-
# file.path(output_fig_wd, 'pcrit-100c') if
# (!dir.exists(output_fig_pcrit_100c_wd)) {
# dir.create(output_fig_pcrit_100c_wd) } ids <- labchart_tidy_fish_100c %>%
# dplyr::distinct(id) %>% pull(id) %>% as.list() # Open a single PDF device
# once pdf(file = file.path(output_fig_pcrit_100c_wd,
# 'combined_pcrit_plots.pdf'), width = 8, height = 6) for (id_i in ids) {
# id_name <- id_i mo2_data <- labchart_tidy_fish_100c %>% dplyr::filter(id ==
# id_i) MR_set <- mo2_data$SMR[1] %>% as.numeric() tryCatch({ # Generate and
# render the plot plot_pcrit( po2 = mo2_data$o2, mo2 = mo2_data$MO2, method =
# 'All', avg_top_n = 1, level = 0.95, iqr = 1.5, NLR_m = 0.065, MR = MR_set,
# mo2_threshold = Inf, return_models = FALSE, showNLRs = FALSE ) # Add a title
# in the top-left corner mtext(text = paste(id_name), side = 3, line = 2, adj =
# 0, # Top margin, aligned to left col = 'blue', font = 2, cex = 1.2) }, error
# = function(e) { message('Skipping channel ', id_name, ' due to error: ',
# conditionMessage(e)) }) } # Close the PDF device *after* the loop dev.off()
Plotting in the html. None of the models appear to estimate a
Pcrit value convincingly.
# ids <- labchart_tidy_fish_100c %>% dplyr::distinct(id) %>% pull(id) %>%
# as.list() for (id_i in ids) { id_name <- id_i mo2_data <-
# labchart_tidy_fish_100c %>% dplyr::filter(id == id_i) MR_set <-
# mo2_data$SMR[1] %>% as.numeric() tryCatch({ # Generate and render the plot
# plot_pcrit( po2 = mo2_data$DO, mo2 = mo2_data$MO2, method = 'All', avg_top_n
# = 1, level = 0.95, iqr = 1.5, NLR_m = 0.065, MR = MR_set, mo2_threshold =
# Inf, return_models = FALSE, showNLRs = FALSE ) # Add a title in the top-left
# corner mtext(text = paste(id_name), side = 3, line = 2, adj = 0, # Top
# margin, aligned to left col = 'blue', font = 2, cex = 1.2) }, error =
# function(e) { message('Skipping channel ', id_name, ' due to error: ',
# conditionMessage(e)) }) }